Exact solution of the Bose-Hubbard model on the Bethe lattice 
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The exact solution of a quantum Bethe lattice model in the thermodynamic limit amounts to 
solve a functional self-consistent equation. In this paper we obtain this equation for the Bose- 
Hubbard model on the Bethe lattice, under two equivalent forms. The first one, based on a coherent 
state path integral, leads in the large connectivity limit to the mean field treatment of Fisher et 
al. [Phys. Rev. B 40, 546 (1989)] at the leading order, and to the bosonic Dynamical Mean Field 
Theory as a first correction, as recently derived by Byczuk and VoUhardt [Phys. Rev. B 77, 235106 
(2008)]. We obtain an alternative form of the equation using the occupation number representation, 
which can be easily solved with an arbitrary numerical precision, for any finite connectivity. We 
thus compute the transition line between the superfluid and Mott insulator phases of the model, 
along with thermodynamic observables and the space and imaginary time dependence of correlation 
functions. The finite connectivity of the Bethe lattice induces a richer physical content with respect 
to its infinitely connected counterpart: a notion of distance between sites of the lattice is preserved, 
and the bosons are still weakly mobile in the Mott insulator phase. The Bethe lattice construction 
can be viewed as an approximation to the finite dimensional version of the model. We show indeed 
a quantitatively reasonable agreement between our predictions and the results of Quantum Monte 
Carlo simulations in two and three dimensions. 

I. INTRODUCTION 

The bosonic version of the Hubbard model was introduced for the first time in the seminal paper [1]. It describes 
a system of bosons hopping between neighboring sites of a lattice, and subjected to a local repulsive interaction 
disfavouring the multiple occupancy of a site. The competition between these two effects leads to a quantum phase 
transition [2] at zero temperature: when the hopping term dominates the groundstate is a Bose-Einstein Condensate 
(BEC) delocalized over the lattice, also known as a superfluid (SF) phase. On the contrary for large local repulsion 
it becomes energetically favorable to form a commensurate state where the average number of bosons per site is fixed 
to an integer value. This incompressible phase is called a Mott insulator (MI). 

The introduction of this model was motivated by experimental results on Helium adsorbed on disordered substrates. 
The recent progresses of the experimental techniques for the manipulation of cold atoms, and in particular the 
possibility of devising optical lattices loaded with cold gases, revivified the interest for the Bose-Hubbard model. 
Following the proposal made in [3] , the experimental observation of this Mott transition was first achieved in [4] . We 
refer the reader to [5] for a review of such experiments bridging a gap between atomic and condensed-matter physics. 

From a more theoretical point of view the Bose-Hubbard model has been studied with various techniques, notably 
mean- field approximations [1, 6], perturbative expansions in the hopping strength [7-9], Random Phase Approxi- 
mations (RPA) [10-14], and Quantum Monte Carlo simulations [15-20]. The mean-field approach of [1] yields a 
qualitatively correct prediction of the phase diagram of the model. However its description of the Mott insulator is 
oversimplified, the hopping being completely neglected in this phase at the mean-field level. As often the assumptions 
underlying the mean-field approximation become valid when the connectivity of each site (the number of its neigh- 
bors) is very large, either of the order of the system size itself (in which case a rigorous treatment of the model is 
possible [21, 22]), or in the limit of large dimensionality. The latter case was recently investigated in [23, 24] where a 
bosonic version of the Dynamic Mean-Field Theory [25] was developed, including first order corrections in the inverse 
of the spatial dimension. 

In this paper we follow a related but slightly different road to improve on the mean-field treatment of [1], by treating 
the Bose-Hubbard model at the level of the Bethe approximation, i.e. by solving it on a Bethe lattice. We show that, 
in the thermodynamic limit, the computation of all observables amounts to solve a single self-consistent functional 
equation. By writing this equation in the coherent state basis one recovers the B-DMFT in the large connectivity 
limit [23, 24]. Unfortunately, this equation cannot be analytically solved for finite connectivity, except in some special 
limits. Hence, we show how to rewrite it in a more convenient way using the occupation number basis; this opens the 
way to its resolution through an efficient numerical algorithm. 

The Bethe approximation is the next step in a hierarchy of approximations for lattice systems known as the Cluster 
Variation Method [26]; it is in fact exact when the model considered is defined on a Bethe lattice, that is a graph which 
is locally a tree (no short loops can be closed on it). Bethe lattices appear in the DMFT derivations [23-25] as an 
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intermediate step, before taking the limit of large connectivity. On the contrary in the present paper the connectivity 
of the Bethe lattice will be kept finite. From the point of view of the imiversality classes of critical exponents Bethe 
lattice models fall in the mean-field category, yet their finite connectivity makes them richer than the (fully-connected) 
mean- field version of [1]. In the latter any site is adjacent to all others, whereas on a Bethe lattice one has a well 
defined notion of distance between two sites, though it does not correspond to the usual Euclidean distance. Moreover 
the MI phase is non-trivial for the Bethe lattice version of the model, at variance with the mean-field picture. Though 
we only consider here the ordered version of the model a motivation for studying the Bethe lattice is the possibility 
that it exhibits a Bose Glass phase [1, 27] in presence of disorder, which cannot be described at the simplest mean- field 
level. This hope is backed up by the relevance of the Bethe lattices for the localization problem [28]. 

The methodology we use has its roots in the extensive research effort which took place over the last decade in the 
community studying statistical mechanics of disordered systems. Under the name of cavity method a set of techniques 
has been developed to solve classical spin models on Bethe lattices (in the sense of sparse random graphs), with 
applications to spin glasses [29] and to random optimization problems arising from theoretical computer science [30] . 
An extensive presentation of this field can be found in the recent book [31]. The cavity method has been recently 
extended from classical to quantum spin models in presence of a transverse field [32-34] . We widen here the range of 
applicability of the method by including lattice boson models. 

The rest of the paper is organized as follows. In Sec. II we define precisely the model studied and the physical 
observables of interest, and we recall the usual mean-field treatment and its qualitative physical predictions. Sec. Ill 
is an overview conveying the main ideas of the cavity method, starting for pedagogical reasons from the ferromagnetic 
Ising model before explaining its extension to quantum problems, and in particular its connection with the B-DMFT 
of [23, 24]. Sec. IV contains the core of our contribution; we first present in Sec. IV A the main equations describing 
the Bose-Hubbard model on the Bethe lattice, and the principles of their numerical resolution. In Sec. IV B we 
present our predictions for several physical observables and compare the phase diagram we obtain with the Quantum 
Monte Carlo results in two [17] and three [18] dimensions, as well as with the B-DMFT prediction of [24] in three 
dimensions. For the convenience of the reader not interested in technical aspects we postpone the detailed derivation 
of the equations to Sec. V. Finally we present our conclusions and perspectives for future work in Sec. VI. 

II. DEFINITIONS AND REMINDER ON THE TRADITIONAL MEAN-FIELD APPROACH 

A. Definition of the model and of the observables 

We shall consider the Bose-Hubbard model for spinless bosons on a graph of N vertices (sites) defined by the 
following Hamiltonian, 

U ^ w 

H = - J ^(atoj + atoi) + ^ ~ = Hk + Hl , (1) 

i=i i=i 

where the first sum nms over a subset of the N{N — l)/2 possible edges (links) between pairs of sites, J is the hopping 
amplitude between neighboring sites, a,; (rcsp. aj) is the boson annihilation (resp. creation) operator on site i, fi is 
the chemical potential fixing the density of particles, and U (;ontrols the strength of the on-site interaction between 
particles. It is convenient to separate the kinetic term Hk (proportional to J) in the Hamiltonian, from the local term 
Hl including the Hubbard interaction and the chemical potential term. Note that the kinetic energy defined in such 
a way will be negative: indeed, the discrete lattice version of the kinetic energy would be given by Hk + JY^i Ciajai 
where Cj is the connectivity of site i. 

Although the method we will discuss in this paper allows in principle to compute very general observables (such as 
multi-point imaginary time correlations), in the following we will be mainly interested in standard observables such 
as the mean density, kinetic and on-site energy, condensate fraction, and Green functions. For the sake of clarity we 
recall now their definitions. The partition function at temperature T is defined by 

Z = Tr[e-f^"] , (2) 

where /3 = 1/T (we set ks = 1), the corresponding free energy is F = —TlogZ; the free-energy per site is / = F/N. 
The thermodynamic average of an operator O is defined as 



(0) = ^Tr[Oe-^^] 



(3) 
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In particular the average number of particles on site i reads (rii) = ^afoiy and the density of particles is p = 

N^^ (rii); the kinetic (resp. local) energy per site is ck = {Hk) /N (rcsp. = (Hl) /N). To define the order 
parameter for BEC, one possibility that is convenient for our purposes is to introduce a small perturbation to the 
Hamiltonian in the form = H + e X^j(ai + a|) and define 

(a) = lim (a), , (4) 

where (•)^ denotes the thermodynamic average in presence of the perturbation. 

The imaginary-time evolution of an operator O is given by 0{t) = e^^Oe~^^; we then define the two-times 
correlation functions as [35] 



Glir) = (a,(r)a|(0)) = ^TT[e 



Z 

The Green function is defined, for —(3/2 < r < /3/2, by 

G\t) = e{T)GUr)+e{-T)G^^{T) ^ (^Ta,{r)al{0)) , (6) 

where T is the usual time-ordering operation, which should not be confused with the temperature. Note that the 
cyclic property of the trace implies that G^(r) = G>(r + hence 

G\t) = e{T)GUr) + 0{-t)GI{(3 + r) , (7) 

and the knowledge of Gl, (r) for < r < /3 is enough to determine the Green function. 

The cavity method allows also the computation of spatial correlation functions, we shall in particular determine 

the one-particle density matrix pij = (ajaj 



B. Review of mean- field results 



The inexistence of an analytical solution of the model for finite dimensional lattices has triggered a large amount 
of research effort on numerical simulations [15-19], pertubative expansions [7, 8], or mean-field treatments [1, 6] of 
the problem. Let us briefly recall the various points of view and methodology that the mean-field approach usually 
encompasses. It can first be seen as an approximation to finite-dimensional models. Maybe the most satisfactory 
way to perform this approximation is by means of a variational method. One introduces a simpler trial (Gutzwiller) 
Hamiltonian where the sites are decoupled. 



N r 



i=l 



U 



(8) 



the V'i being here free parameters. Note that i7o breaks the particle conservation symmetry. The free-energy at 
inverse temperature (3 of the original model, F(i3), can be upper-bounded as F < Fo -I- (i?)o — {Ho)o, where Fq{(3) 
is the free-energy of the trial Hamiltonian and {■)o denotes a thermal average with respect to the trial Hamiltonian. 
The best variational description is thus obtained by minimizing the upper bound with respect to the parameters ipi. 
Assuming that all sites of the graph have the same number c of neighbors (c = 2d for a d-dimensional hypercubic 
lattice), the bound is minimized by taking the same value ip on all sites, which can be chosen real without loss of 
generality. The free-energy per site can then be upper-bounded as /(/3) < /var(/3), with 



/var(/3) = inf 



1 



1 



■ InTr 



(9) 



where the trace is over the Hilbert space of a single site. The physical properties of this variational free-energy are 
well-known [1]: at all temperatures and chemical potentials there is a transition encountered upon increasing the 
hopping intensity J from an "insulating" phase, characterized by (a) = 0, to a "superfluid" phase with (a) ^ 0. At 
zero temperature this transition line draws a series of lobes in the (J/U, ii/U) plane, inside each lobe the particle 
density being constrained to a given integer. These Mott Insulator phases are thus incompressible. There exists 
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alternative ways to obtain the estimation (9) for the free-energy per site of the Bose-Hubbard Hamihonian. The 
usual mean-field approximation consists indeed in replacing ata^ with {al)aj +al{aj) — (ct|)(aj) in the hopping term 
of the original Hamiltonian, neglecting correlations between neighboring sites. This transforms the Hamiltonian into 
the site-decoupled form (8), with tpi given by a sum of (aj) over the neighbors j of i. These expectation values 
are then computed with respect to the decoupled Hamiltonian, the self-consistency equations leading finally to the 
same expression (9) of the free-energy per site as the variational approach; the latter has however the advantage of 
being better controlled, in the sense that it provides a rigorous bound on the true free-energy of the system. Finally, 
another reasoning yielding this mean-field result consists in devising a model which has exactly (9) as its free-energy 
per site in the thermodynamic limit, instead of taking it as an approximation for the finite-dimensional case. As 
could be expected this model corresponds to the fully-connected version of the Hamiltonian (1), with the sum in the 
hopping term running over all possible pairs of sites, with a coupling constant J inversely proportional to the size N 
of the system in order to have an extensive thermodynamic limit. It has been shown rigorously in [21, 22] that in the 
thermodynamic limit the free-energy of this fully-connected model converges indeed to (9). 

The above described mean-field treatment has limitations both of a quantitative nature (the approximation cannot 
be expected to be very precise for small dimensions) and of a qualitative one. In particular the MI phase is rather 
trivial; as the infimum in (9) is reached in ip = 0, the hopping of particles is completely suppressed in this phase. This 
drawback is particularly severe in the case of the disordered Bose-Hubbard model, for which it forbids the description 
of the Bose Glass phase predicted in [1] . In order to cure this defect of the mean-field treatment one has to account 
in some way for the correlations between neighboring sites. Thinking in terms of classical spin models, the mean-field 
approximation is the lowest level of a hierarchy of descriptions (known as Cluster Variation Method [26], or Kikuchi 
approximations [36]) which treats exactly larger and larger regions of the interaction graph. In this work we shall 
deal with the Bose-Hubbard model at the next level of the hierarchy, known as the Bethe approximation, in which 
correlations between nearest neighbors arc explicitly taken into account. In the same way as the simplest mean-field 
approximation was exact for the fully-connected model, the Bethe approximation is exact for a family of graphs, 
known as Bethe lattices. In these graphs each site interacts with a finite number c of neighbors, say 2d in order to 
mimic a d-dimcnsional hypercubic lattice, but the short loops present in the latter are discarded: Bethe lattices have 
a local tree structure, sec Fig. 1 for a picture of the local appearance of a Bethe lattice of connectivity c — 4. For 
at least two reasons it is however better not to picture a Bethe lattice as a finite regular tree (usually called Cayley 
tree in this context). Indeed a regular tree is strongly inhomogenous, a finite fraction of its "volume" being very 
close to its "surface" , and only the center of the Cayley tree has the bulk properties one is interested in. Moreover 
in the case of frustrated models the boundary conditions imposed on the leaves of the Cayley tree have to be chosen 
with particular care. For these two reasons it has become customary in the community of statistical mechanics of 
disordered systems, following [29] , to define a Bethe lattice as a random c- regular graph [37] , that is a graph chosen 
uniformly at random from the set of the graphs on N vertices where all sites have precisely c neighbors. These Bethe 
lattices are locally tree-like, their loops have typically a length diverging with the size N of the system, yet they 
do not have any boundary, all sites playing the same role (in the same way as periodic boundary conditions impose 
translation invariancc on a finite cubic lattice). The absence of an underlying finite-dimensional structure gives them 
a mean-field character, but their finite connectivity leads to a richer content than the fully-connected models. The 
free-energy of Bethe lattice models is self-averaging with respect to their random character in the thermodynamic 
limit. In other words for large enough N a single sample is a good representative of the ensemble average. The 
so-called cavity method has been developed to treat classical spin models defined on such random graphs, and found 
important applications for optimization problems of theoretical computer science [30] . 



For the sake of pedagogy we shall begin our presentation by the cavity method treatment of the ferromagnetic Ising 
model on the Bethe lattice. We consider the following energy function. 



where the classical degrees of freedom ai can take the values ±1, the first term describing ferromagnetic (J > 0) 
interactions between neighbors on a graph of N vertices, the second one corresponding to an uniform magnetic field. 
We shall repeatedly use in the following the notation di for the set of vertices adjacent to a given vertex i, i.e. for the 
sites which interact with i, and di\j for those vertices around i distinct from j. The Gibbs-Boltzmann probability 



III. OVERVIEW OF THE QUANTUM CAVITY METHOD 



A. The cavity method for ferromagnetic Ising models 




(10) 
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FIG. 1: A portion of a Bethe lattice of connectivity c = 4, approximation of a square lattice. 



measure is 



r?(ai,...,c7;v) = |e-'5^(-^ , (11) 

with the partition function Z ensuring its normahzation. The goal is to compute the free-energy per site / = 
— {\nZ)/{N(3) and the local magnetizations ((Ji), the brackets denoting an average with respect to the law (11). 

Let us first consider the case where the interaction graph is a finite tree. In this case the computation can be 
organized in a very simple way, taking benefit of the natural recursive structure of a tree. We define the quantity 
Zi^j{ai), for two adjacent sites i and j, as the partial partition function for the subtree rooted at i, excluding the 
branch directed towards j, with a fixed value of the spin variable on the site i. We also introduce Zi{ai), the partition 
function of the whole tree with a fixed value of o-j. These quantities can be computed according to the following 
recursion rules, see Fig. 2 for an example, 

kedi\j \ o-fc / j^di \ (Tj J 

It will be useful for the following discussion to rewrite these equations in terms of normalized quantities which 
can be interpreted as probability laws for the random variable ai, namely r]i^j{ai) = Zi^j{ai)/^^, Zi^j{a') and 
Vi{(^i) = / J2a' ^he recursion equations read in these notations 

^'^^ kGdi\j \ / jedi \<T, J 

where Zi^j and Zi are normalization constants. On a given tree the recursion equations on the rji^j for all directed 
edges of the graph have a single solution, easily found by propagating the recursion from the leaves of the graph. 
Moreover the quantity r]i{ai) is exactly the marginal probability law of the Gibbs-Boltzmann distribution (11), hence 
the local magnetizations can be computed as {ai) = r]i{a)a. 

The reasoning above was made under the assumption that the interaction graph was a tree. One can however 
look for a solution of the recursion equations (13) on any graph, even in the presence of loops. This approach is 
known as Belief Propagation in inference problems [38], and corresponds to the Bethe approximation of statistical 
mechanics [39] . The cavity method is a set of prescriptions to use these local recursion equations to predict the behavior 
of models defined on sparse random graphs, that are locally tree-like (the typical length of the loops diverging in the 
thermodynamic limit). In its simplest version (known as the replica symmetric one) one makes the assumption of 
the existence of a single pure state which implies that the effect of the loops does not spoil the local recursions (13). 
Their presence simply provides self-consistent boundary conditions and avoids the ill-defined behaviour due to the 
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FIG. 2: Example of an Ising tree model on 7 vertices. The definition of Z'z^x and its recursive computation reads here: 



dominant surface effects of trees. In the case of unfrustrated ferromagnetic models it has been shown rigorously [40] 
that this assumption is correct, the predictions of the cavity method being exact in the thermodynamic limit, both 
for the local magnetizations and for the free-energy per site. Let us be more explicit for the case of the Bethe lattice, 
where all vertices have the same connectivity c. The probability laws rji^j are then the same on all edges; denoting 
r]ca,v their common value, one finds the self-consistent equation 

»7cav(^) = — e'^'*^ ^-v(ai) . . . %av(c7c-l)e^^^("^ + -+"=-^) , (14) 

with ^cav a normalization constant. A pictorial representation of this equation can be found in Fig. 3. The local 
magnetization is then computed as 

(a) = ^ i^{u)a , r,(a) = ie'^"- ^ r,eav(ai) . . . r;eav(ac)e^^-('^^+-+'^=) , (15) 

cr (Ti,...,(Tc 

including the c neighbors of a central site as represented in Fig. 4. The term cavity comes from the fact that in 
Eq. (14) one site has been removed from the neighborhood of the considered vertex. As the Ising variable a can only 
take two values, each of the probability distributions r/cav and rj can be parametrized by a single real number, a cavity 
(resp. effective) magnetic field, 

??cav(cr) = — — , r?((T) = , , , (16) 

2 C0Sh(/j/lcav ) 2 COSh(/j/leflF ) 

solutions of 

c — 1 c 
hca,v = h-\ — arctanh[tanh(/3J) tanh(/3/icav)] , /leff = + "^arctanh[tanh(/3J) tanh(/3/icav)] • (17) 

P P 

thus making the resolution of (14), (15) extremely simple. In particular at zero external field /i = 0, one finds that a 
phase transition occurs at /3 = /3c; separating a high temperature paramagnetic phase (/icav = ^eff = 0) from a low 
temperature ferromagnetic phase (/icavj/ieflf 0)- The critical temperature is easily obtained linearizing the equation 
for /icav around /icav = and is the solution of (c — 1) tanh(/3c J) = 1- 

It is worth noting that one can compute explicitly the spin-spin correlation function, which is given in the param- 
agnetic phase by Cij = {criaj) = [tanh(/?J)]''(*'-'), where d{i,j) is the distance between sites i and j, defined as the 
length of the shortest path connecting sites i and j. The associated correlation length is ^ = — log[tanh(/3J)] which 
is finite at the transition point /3 = (3c- Nevertheless, the associated magnetic susceptibility % = N~^'Y^j^c\ {ai) /dh 
diverges: one can show using the fluctuation-dissipation theorem that x = J2ii ^ij = P Yl'dLo-^dCd, where Nd 

is the number of points at distance d from a given reference point and scales as (c — 1)'' for large d. Therefore, if Cj. 
decays slower than (c— l)^'', the corresponding susceptibility is divergent; this is indeed consistent with the equation 
for /3c obtained above. Hence phase transitions on Bethe lattices are always associated to diverging susceptibilities 
and finite correlations lengths (see [41] for a discussion of this fact in the context of quantum spin models). For finite 
dimensional lattices, Md grows as a power of d, and one needs a diverging correlation length to obtain a diverging 
susceptibility. 
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FIG. 3: Pictorial representation of Eq. (14). The bubbles on the first panel represent subtrees of the graph; their effect on the 
spins ai, . . . ac-i is summarized by r?cav, represented as a bold arrow in the second panel. Tracing over these c — 1 spins leads 
to the third panel. 




FIG. 4: Illustration of Eq. (15); the local magnetization of one site is computed by taking into account all the c neighbors. 



B. Suzuki- Trotter formalism and quantum cavity method 



We have just seen how the cavity method deals with classical spin models on locally tree-like graphs. The extension 
of the method to quantum models is based on the Suzuki- Trotter formula. Generically speaking the resolution of 
a quantum model amounts to compute its partition function Z = Ti [e~^^'\ , where the Hamiltonian H is now an 
operator, which can usually be split as H = Hi + H2 with Hi 2 two operators which do not commute. For instance 
in the case of quantum spin 1/2 models one can have an interaction term Hi involving only the z component of 
the spins, while H2 is a transverse field acting in a perpendicular direction. In the case of the Bose-Hubbard model 
this decomposition splits the Hamiltonian between the hopping term Hk and the on-site potential energy Hl- The 
non-commutativity of Hi ^2 can be cured with the application of the Suzuki- Trotter formula [42], 

Z = lim Tr 

The computation then proceeds with the insertion of Ns representations of the identity between the Ns elements of 
this product. It is convenient to express the identity operator in a basis where one of the two operators is easily 
evaluated. In the case of quantum spin 1 /2 models one can for instance use the eigenvectors of the Pauli matrix in the 
z direction; the quantum partition function is then turned into the partition function of a classical Ising spin model, 
where each quantum spin is replaced by a set of Ng classical spins, indexed by their position on the original graph and 



-Hi 



(18) 
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a supplementary "discrete imaginary time" coordinate (which becomes a continuous parameter in the iVg ^ oo hmit). 
The important point is that the spatial structure of the graph of interactions is preserved in this construction, at the 
price of the introduction of imaginary time-dependent classical degrees of freedom. In particular, if the interactions 
of the quantum model are defined on a tree-like graph, the cavity method still applies to this extended classical 
model. This line of thought was first followed for quantum spin models in [32] (see also [33]) for a finite number Ng of 
Suzuki- Trotter slices, the continuous imaginary time limit was then taken in [34]. In this paper we shall extend this 
method to deal with lattice bosons models. 

In this case the decomposition of the identity operator in the Suzuki- Trotter can be expressed using either coherent 
states or occupation numbers. The latter has the advantage of being discrete, and we shall use it in the rest of the 
paper. For the sake of clarity and in order to make contact with the recently proposed B-DMFT [23, 24] we discuss 
first the application of the cavity method within the coherent states representation in the rest of this subsection. 
Inserting such a decomposition of the identity for each of the sites at each of the Suzuki- Trotter slices leads, 
in the continuous time limit Ng —^ oo, to the coherent state path integral expression of the partition function of the 
Bose-Hubbard model [35]: 



N 

' /n 

i=l 



DttiDtti e ^ 



dr 



N 



u 



ai{T){dr - M)ai(T") + -wiai{T-)ai{T)f ) - 



{i,3) 



(19) 



(20) 



Here and in the following we use bold symbols to denote imaginary-time dependent quantities. ai{T) and ai{T) are 
two (formally conjugate) complex numbers indexing the coherent state at site i and imaginary time r, with DaiDdi 
the path-integral measure of this site. Following the same steps as in the study of the Ising model, the cavity method 
for a Bethe lattice of connectivity c leads to the following self-consistency equation: 



,(a, a) = w{a,a) / T] DaiDaiT]cav{ai,di) exp J / dr 

Zcajv J .I \ Jo 



air) ai(r) + a{T) ^ a,(T) 



(21) 



with the on-site weight of the path (a, a) given by 



'w{a, a) = exp 



-£dr (a(r)(a. -M)a(T) + ^{a{T)a{T)f'^ 



(22) 



This self-consistent equation on Tjcav is formally similar to the corresponding Eq. (14) of the Ising model. It is 
however much more complicated: the Ising degree of freedom a could only take two distinct values, whereas (a, a) 
belongs to a space of hmctions of the imaginary time t. Therefore r;cav(fl, o) is a functional measure whose repre- 
sentation is much more difficult; a complete parametrization requires the knowledge of all n, m-points correlations 
{a{ti) ■ ■ ■ a{tn)a^{si) ■ ■ ■ a^Sm))- On the other hand, this is one of the most interesting features of the quantum cavity 
method: on-site quantum fluctuations are fully kept into account, without any truncation of higher order correlations. 

Unfortunately, an exact solution of (21) can be easily obtained only in the case of free bosons {U = 0). rjcav acquires 
in this case a Gaussian form, with averages and two-point functions which can be computed exactly and reproduce 
the results obtained by direct diagonalization of the adjacency matrix of the Bethe lattice [43] . 



C. Large connectivity limit and the connection with B-DMFT 



In the interacting case (J7 > 0) a solution of (21) can be looked for in the limit of large connectivity, and this is 
precisely the road followed by the B-DMFT studies [23, 24] . For completeness we shall explain in this subsection how 
to recover the B-DMFT formalism from Eq. (21), before turning in the next section to the occupation number basis 
which will allow to solve the model for any connectivity. 

To state the large c expansion let us rewrite (21) in a more convenient way, using the lighter Nambu notation 
^^{t) = {d{T),a{T)), and consequently = (a, a). Then we can rewrite (21) as 



,(*) = J_ w{^) e('^-i)r(-^*) , r(*) = log I / £)* 77cav(*) exp dr ^^{t)^t) 

•2'cav \J Jo 



(23) 
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r(*) is the generating functional of the connected correlation functions of a and a. It can be expanded as 

r(*)=/ dr + / dTdr'$t(-r)Ge(T-T')4'(r') +0(4'^) , (24) 

Jo Jo 

where the averages (•) arc with respect to r^cav, we have used the cyclic invariance in imaginary time, and Gc{t — t') = 
(\E'(t)5'^(t')) — {^) (^^) is the connected part of the two point correlator of ^ . 

In the large connectivity limit, the supcrfluid- insulator transition happens for a critical value oi J = / with a 
finite . This can be argued by looking at Eq. (9), where it is clear that the dependence on hopping and connectivity 
is only through J = Jc, which is the real control parameter. For large c and J = J /c, the cumulant expansion of 
(c — l)r(J'^/c) thus becomes a systematic expansion in powers of 1/c. By keeping only a finite number of terms in 
the cumulant expansion, wc obtain an expression of ?]cav('5') that is not Gaussian (because of the U term in w(^')); 
still, we can obtain closed equations for the cumulants by computing them self-consistently as averages over the 
non-Gaussian r/cav 

The leading order in c gives, assuming without loss of generality that the average value of the order parameter is 

real, (*t) = (^,^), 

^^^^(^) = ^^(^)e/o''dr^(*t)*(.) ^ i,= — j i)a£>aa(0)e-/o''<i^(»(^'(^^-'')»(^)+*(»(^)»(^))'-'^^(«M+"M)) . 

•^cav -^^cav J 

(25) 

This last equation can be rewritten in the operator representation, which gives back the equation for ip corresponding 
to the minimization of the variational free-energy (9), up to a multiplicative constant in the definition of tf). Note 
that a generalization of the discussion above and of Eq. (25) to the disordered case leads to the stochastic mean field 
theory devised in [44]. 

We will now show that the next-to-leading order in the cumulant expansion of (24) gives the B-DMFT equations 
recently derived in [23, 24]. Note that the truncation at this two-point level was also used in the context of spin 
models in [32]. Plugging the expansion (24) in (23), we obtain to order 1/c: 

r/cav(*) = exp [ - S\oc] , 

•^cav 

Sxoc= dTdT'^\T)g-\T-T')^{T')+ [\t \^i^\T)^{T)f - {^^)^{t)\ , (26) 

Jo Jo 1° c _ 

G-Hr -r') = l - ^ ^) Sir - r') - ^S.(r - .') . 

Then, (^'^) and Gc(t — r') = {^{t)^"^ {t')) — {^) (^'^) have to be computed self-consistently as averages with the 
local action 5ioc- This set of equations correspond exactly to the B-DMFT of [23, 24] for the special case of a Bethe 
lattice. 

Away from these two limits {U = Q and c oo) it seems difficult to obtain a solution of the cavity equation (21) 
as written in the coherent state basis. As a consequence we shall turn in the following to the representation number 
basis to apply the Suzuki- Trotter formula and thus obtain a more tractable equation for all values of U and c. 



IV. THE QUANTUM CAVITY METHOD IN THE OCCUPATION NUMBER BASIS 

A. The equations and the procedure for their numerical resolution 

The insertion of a decomposition of the identity expressed in the occupation number basis in the Suzuki- Trotter 
formula (18) leads to an expression of the partition function of the Bose-Hubbard model as a sum over occupation 
number trajectories in imaginary time, {^^(t)}. These trajectories are defined on an imaginary time interval of 
length /3, with the periodicity condition ni(0) = ni{j3). The weight (action) of these trajectories has two origins: the 
local part of the Hamiltonian (1) yields a contribution of the form exp[— j dTV{ni{T))\ for each of the sites, where 
V{n) = Un{n— 1)/2 — fin is the local energy term in the Hamiltonian. In addition the hopping term of the Hamiltonian 
imposes constraints between the occupation number trajectories: each time ruir) is raised (resp. decreased) by 1, 
the occupation number nj{T) of one of the neighbors j € di must decrease (resp. increase) of 1, meaning one particle 
has jumped from j to i (resp. from i to j). Moreover each hopping event multiplies the weight of the trajectory 
{ni{T)} by J and by a coefficient depending on the instantaneous occupation numbers of the sites involved in the 
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hopping. An explicit derivation of this representation shall be given in Sec. V. Here we wish to present the results 
of the computation in a lighter way for the ease of the reader not interested in the technical details (we follow the 
methodology developed for quantum spin models in [34]). 

The recursion equations of the Bose-Hubbard model defined on a tree (or on a locally tree-like graph with the 
assumptions of the cavity method) can be expressed in terms of the probability distribution ?7cav(y) of the hopping 
trajectories y defined on the edges of the graph. This quantity y encodes the imaginary times at which a particle has 
crossed the edge, and the directions of the jumps; examples are given in Fig. 5. In the case of a regular Bethe lattice 
of connectivity c one obtains the following self- consistent equation on ?7cav(y), 

»7cav(y) = — ^"link(y) X! '7cav(yi)---??cav(yc-l)"'iter(y,2/l,---,yc-l) • (27) 

^^^^ 

For simplicity we denote here with a sum symbol what is actually an integral over continuous degrees of freedom y. 
The explicit expressions for the weights w shall be given in Sec. V, see Eqs. (46), (47). Let us emphasize the formal 
similarity with the Ising equivalent Eq. (14), and of course the greater complexity of the basic degree of freedom 
in the Bose-Hubbard case, the hopping trajectory y assuming values in a much larger space than the Ising variable 
(J e {+1,-1}. We already mentioned this problem while discussing the related equation (21) in the coherent state 
basis. Compared to this latter case we are however facing now an easier problem: the number of hopping events on 
a given edge is a (random) number which remains finite as long as the temperature is positive (it is actually of order 
(3 J). A single hopping trajectory y can thus be encoded as an integer p, i.e. the number of particle jumps occuring 
on this edge during the imaginary time interval [0, p reals precising the imaginary times where these jumps occur, 
and p binary variables giving the direction of the jumps, see Fig. 5. In contrast the coherent state trajectories were 
those of two continuous functions with a priori no compact representation. This remark opens the way to an efficient 
numerical method for the resolution of (27) [60]. We can indeed follow the population dynamics strategy [28, 29]. 
The idea of this method is to represent numerically the probability distribution r]cav{y) as a (weighted) sample of a 
large number ^traj of hopping trajectories, namely 

??cav(y) ^ giSiy -Vi) ^ (28) 

i=l 

where the A/traj weights of the trajectories are normalized according to 

Sampling an element y from the probability distribution ?7cav corresponds in this representation to extract an integer 
i € [1, A/traj] with probability gi, and setting y = y^. This representation of Tjcav is an approximation, which yields 
better and better numerical accuracy when A/traj grows. The determination of a sample of weights gi and trajectories 
y^ which turns the representation (28) into a good approximation of the solution of (27) can be performed iteratively. 
To explain this point let us first rewrite the self-consistent equation (27) as 

?7cav(y) = X] ^cav(yi) • • ■ ??cav(yc-i)-P(y|yi, • • • ,yc-i) ^^^^' ' " " , (30) 

1/lviS/c-l 

where we have defined 

D/ I ^ ■'^link(y)witor(y,yi,.-.,yc-l) ^, \ ^ / ^ / \ lr,^\ 

P{y\vi,- ■■,Vc-i) = 57 ^ ' ^(yi,--->yc-i) = >^wiink(y)witer(y,yi,..-,yc-i) • (3i) 

^(yi,---,yc-i) ^ 



y 



Constructed in this way -P(y|yi, • • • ,yc-i) is a conditional probability distribution over y. Given the values of the 
hopping trajectories y^, . . . ,yc_i it is actually possible to perform an exact sampling from P{y\yi, . . . ,yc-i) ^"^^ to 
compute the normalization constant Z{yi, . . . ,y^_i). In fact this reduces to a relatively simple single site problem: 
one has to construct the occupation trajectory n^r) of a site, and the associated hopping trajectory y towards one 
of its neighbors, given the hopping trajectories yi, . . . ,yc-i on the other adjacent edges (see Fig. 5 and recall the 
pictorial representation given for the Ising case in Fig. 3). The c — 1 hopping trajectories on the upper edges impose 
that n(T) changes by ±1 at the times particles arrive or depart from the considered central site. Between these times 
we shall show in Sec. V that the single site problem is described by an effective Hamiltonian V{n) — J{a + a^); each 
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y2 = (0) 
y3 = (lT^) 

y=(Xrr, irr) 
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J 



1 



FIG. 5: Example of the sampling process from the law P{y\yi, 
hopping trajectories (j/j, i/j; J/a) impose jumps in r{" < < 



• • • , Vc-i) defined in Eq. (31) for c = 4. The c — 1 incoming 
< rl", a particle arriving on the central site at the first three 
imaginary times, leaving at the fourth one, corresponding to the symbols in the upper part of the rightmost figure. A trajectory 
n(T) respecting these constraints is shown; it contains two other jumps at rf"* < r2"', associated to hopping events in y. The 
symbols on the bottom are reversed: particles leaving the central site will arrive on the next level of the recursive equation. 



change in the value of ^(t) provoked by the creation/annihilation operators in this effective Hamiltonian is associated 
to an hopping event in the new (downwards) trajectory y. An example of this construction is given in Fig. 5. The 
computation of Z{y-^^, . . . ,yc-i) can be performed recognizing it as the partition function of the single-site effective 
Hamiltonian. 

Suppose now a representation of r^cav is available under the form (28): one can phig it into the right-hand-side 
of Eq. (30) and construct a new set of trajectories and weights corresponding to its left-hand-side, let us call them 
{yjj9j}flT'- Independently and identically for each j G [l,M;raj]j this can be done by: 

• drawing c — 1 integers ii, . . . , ic-i € [1, A/traj], each of them independently with probability gi 

• drawing y from P{y\y^^ , . . . , y,^_ J 

• setting y'^ = y, g'^ = Z{y^^ , . . . , _ J 

The new weights g'^ are then normalized to fulfill Eq. (29). This process can be iterated, plugging in the r.h.s. of 
Eq. (30) the newly obtained representation {y^ , 51^ }, and so on and so forth. After a certain number of steps, starting 
from an arbitrary initial condition [61], a stationary solution is reached. 

The physical obsorvables can be computed from this representation of r?cav Let us first explain the determination of 
the average occupation number of one site. Following our convention we denote n the occupation number imaginary- 
time trajectory ^(t). Its probability distribution rj{n) is expressed by considering the complete environment of a 
vertex with its c neighbors, as was done for the local magnetization of the Ising model in Eq. (15) and schematized 
in Fig. 4, 

77(n) = — ^ ?7cav(yi)...??cav(yc)w'site(n,yi,---,yc) • (32) 

The explicit expression of the weight Wsite shall be given in Sec. V. It is however intuitive that it will contain a 
factor cxp[— J dTT^(n(r))] arising from the local part of the Hamiltonian, and requires a consistency between n and 
[y^, . . . ,y^. In fact all the discontinuities in the occupation trajectory n(T) are fixed by the hopping events on the 
c neighboring edges. If the number of hoppings towards the considered site equals the number of jumps outside, 
then the trajectory r?.(r) is fixed up to a global shift ^(t) "> ^(t) -|- m with m independent of time. Whenever 
there is an unbalanced number of hopping towards/outside the vertex no periodic trajectory ^(t) can be constructed, 
Wsito = for such a configuration of (yi, • • • ,yc)- The average occupation number of a site is easily obtained from 
this probability ry(n), 

E Vce^Avi) ■ • •%av(yc)Ew'site(n,yi, • • ■,yc)n{0) 

(a'^a) = r](n)n(0) = '"'.^ 7 — r 7 — 7 r — 

n ^ ??cav(2/l)---%av(yc)Ew'site(Tl,yi,---,yc) 



(33) 
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FIG. 6: (Main panel) (n) and (a) as functions of ^/U for J/U = 0.02, at f3J — 1 and TVcut = 6 for Bethe lattices of connectivity 
c = 4. Full lines are obtained by joining points measured every Afi/U ~ 0.05, except close to the critical points where 
A^/U ~ 0.001. (Inset) Blow-up of the region close to the first SF-MI transition that happens at /Xc = 0.0648 [/. Here (o) is 
plotted as a function of 5 = (/ic — fJ-)/U (circles). The full line is a fit to (a) ~ A5°-^ with A = 2.9475. 

where in the last expression we have explicited the normahzation constant ^sitc. Note that the time r = where 
n(T) is evaluated is arbitrary because of the cyclic invariance along the imaginary time axis, hence n(0) can be 

equivalently replaced by the temporal average ^ Jq dTn(T). Given a representation of yycav as a weighted sample (28) 
the numerical determination of {^a^aj is straightforward. Both the numerator and the denominator of (33) have the 
form of the average of a function u(yi, . . . with the y's independently drawn from their distribution ?7cav As 
already explained drawing from r/cav a trajectory y amounts to draw an index i S [l,A/traj] with probability gi and 
picking the element of the population. In formula this method of sampling leads to 



where {ij,/} are integers drawn independently with the probability gi, with j G [1, A/trics] and Z £ [Ij c]. The numerical 
accuracy of the sampling is increased by taking a large value of A/trics. Moreover one can interleave average steps with 
iteration steps, determining a new set of weights gi and trajectories y^ and recomputing independent averages on this 
new representation of rycav • 

We shall show in Sec. V that all the other observables defined in Sec. II (order parameter (a), free-energy, ki- 
netic/potential energy and Green function) can also be expressed as averages of the form (34), and thus can be 
efficiently computed by the method developed here. 

B. Results 

In the following we present the results obtained for the Bose-Hubbard model on the Bethe lattice solved with the 
cavity method. We have shown in the previous sections that in the thermodynamic limit the exact resolution of the 
model amounts to solve Eq. (27) and we explained how this can be done with arbitrary numerical precision using 
the population dynamics algorithm described in subsection IV A, encoding rjcaviv) as a population of a large number 
A/traj of trajectories. We typically used Atraj = 32768 and checked that the results were unchanged when using larger 
Atraj- We recall that in the thermodynamic limit the local observables (say, {ui}) are independent on the site i. In 
the following we focus on Bethe lattices of connectivity c = 4 and 6, which mimic the connectivity of the two and 
three-dimensional square and cubic lattices, respectively. 

Note that in order to implement the summations over n involved in the computation of the various weights appearing 
in Eqs. (27) and (32), we have to introduce a cutoff on the possible values of n, < n < A^cut- The technical details 
will be discussed in section V. Since all the results presented below corresponds to regions of the parameter space 
such that < (n) < 3, we work with a cut-off on the occupation number in the interval 4 < A^cut < 6. We have 
checked for various values of the parameters that all the results are stable with respect to changes of the cutoff (in 
some cases also up to iVcut > 6). 
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FIG. 7: (Main panel) Total energy, e = {H)/N, as a function of n/U for J/U = 0.02, at I3J = 1 and A^cut = 6 for Bethe lattices 
of connectivity c = 4. (Inset) Kinetic energy ex as a function of the chemical potential for the same values of the parameters. 



In this section we also compare our results with other analytic approaches, such as the variational mean- field 
treatment described in Sec. II B and the bosonic Dynamical Mean Field Theory (B-DMFT)[62] discussed in Sec. IIIC, 
that correspond to the large c limit of our method [23, 24]. We also compare them to Quantum Monte Carlo simulations 
on 2d square lattice [17] and 3d cubic lattice [18]. 

1. Thermodynamic observables 

Once the numerical iterative procedure has converged to a fixed point solution, one can compute the thermodynamic 
observables. Fig. 6 shows the behavior of (n) and (a), as a function of the chemical potential, ^/U for J/U = 0.02, 
and PJ — 1 for a Bcthc lattice of connectivity c = 4. We have also computed (n) and (a) at lower values of the 
temperature {(3J = 2, 4, and 6), but we find that the results are unchanged within our numerical accuracy, except 
very close to the lowest-/i superfluid-insulator transitions, where a little temperature effect is detected[63] between 
(3J = 1 and (3 J = 2 (while no change is detectable above /3J = 2). Hence, we can safely assume that for (3J > 2 we 
are in the zero temperature regime. 

The plot of Fig. 6 clearly shows a sequence of transitions from superfluid phases (SF) to Mott insulators (MI). 
The superfluid phase is characterized by a finite value of (a), which corresponds to a finite condensate fraction 
/bec = («)^/ {n-), while in the Mott insulator (a) = 0. In the limit of zero temperature in the MI regions the average 

number of bosons per site, (n) is fixed to integer values, (n) = 1, 2, 3, As a result, the compressibility of the 

system, x = d{n)/dii vanishes identically. At small but finite temperature, we expect the compressibility to be 
exponentially small in the temperature, and indeed in Fig. 6 wc observe that the compressibility at commensurate 
density is practically zero even at [3 J = 1. 

As shown in the inset of Fig. 6, the critical exponent for (a) ~ (/^c — m)*^ tlie transition from the SF phase 
to the MI is consistent with the mean-field value e = 1/2. As in the classical case, the critical exponents found 
on the Bethe lattice are equal to the mean-field ones (sec [34] for a detailed discussion in the case of the quantum 
spin-1/2 ferromagnet on the Bethe lattice). Therefore, we also expect the other critical exponents to coincide with 
their mean-field values. 

In Fig. 7 we plot the total energy, e = {H)/N, as a fimction of the chemical potential for the same values of the 
parameters as in Fig. 6. We have also computed the free energy of the system, /, but we do not show it since at these 
low temperatures the difference between the energy and the free energy is very small, and practically undetectable 
within our numerical precision, which confirms that we are effectively in the zero-temperature regime and the system 
can be considered to be in its ground state. In the inset of Fig. 7 we show the behavior of the kinetic energy, ck as 
a function of n/U. The kinetic energy is, of course, lower in the SF phase than in the MI. Moreover we find that ck 
does not vanish in the MI phase, differently from the traditional mean-field approach, where hopping is completely 
suppressed in this phase. Note that both the kinetic energy and the potential energy show a discontinuity in the first 
derivative at the transition between the MI and the SF phase. On the contrary, the first derivative of the total energy 
(and of the free energy as well) is continuous at the transition; only its second derivative shows a discontinuity. 
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J/U J/U 

FIG. 8: "Zero temperature" phase diagram of the Bose-Hubbard model on the Bethe lattice in the [J /U, ji/U) plane at 
connectivity c = 4 (left panel) and c = 6 (right panel). The results obtained within the cavity method (black lines and filled 
circles) are compared to Monte Carlo simulation (blue lines and open diamonds) on the square (data from [17]) and the cubic 
(data from [18]) lattice, and to the prediction of the mean-field approach (red dashed curves) and of bosonic DMFT (brown 
full curve in right panel, data from [24]). 



2. Phase diagram 

In Fig. 8 we present the "zero temperature" (recall that our method only works at finite temperature, hence we 
use quotes to remind that an extrapolation to T = is needed) phase diagram of the Bose-Hubbard model on the 
Bethe lattice of connectivity c = 4 (left panel) and c = 6 (right panel) in the ( J/C/, m/C/) plane. The phase boundaries 
have been obtained by computing (a) as a function of at constant J/U. In this case we used 4 < iVcut ^ 6 
and 4 < /3J < 6 (lower values of (5 J and higher values of iVcut have been used at higher ^). As already discussed, 
at these values of /3 we are effectively in the zero temperature regime. The lobelike shape of the phase boundary 
between MI and SF phases is qualitatively similar to the one found in the mean field approach. In Fig. 8 we also 
compare the results found with the quantum cavity method with the outcomes of Quantum Monte Carlo simulations 
(on the square lattice for c = 4 [17] and on the cubic lattice for c = 6 [18]), and with the results of other analytical 
approaches. This comparison shows that the cavity method does a fairly good job in locating the lobelike contours 
between the MI and the SF phase and performs much better than the mean-field approach. It also performs slightly 
better than the Bosonic DMFT [24] for c = 6, although the difference between the cavity method and the B-DMFT 
is expected to become small as the connectivity is increased. It would be interesting in this respect to compare the 
cavity method and B-DMFT for c = 4. 

In addition to the low temperature phase diagram we have computed the transition temperature line from the 
"insulating" phase (which we defined by (a) =0, but has still a finite conductivity at finite temperature) to the 
superfluid one (defined by (a) ^ 0) at unit filling ((a'I'a) = 1), for Bethe lattices of connectivity c = 6. This result 
is plotted in Fig. 9 and compared to the mean-field prediction and to the three-dimensional Monte-Carlo simulations 
of [18]. 

3. Green functions and particle occupation statistics 

As far as the results shown up to now are concerned, there are no qualitative differences between the Bethe lattice 
and the fully connected models (or equivalently the mean field approximation). We shall now present the results for 
the one-particle on-site imaginary time Green function, which exhibits a richer behavior with respect to the mean-field 
description. This quantity, defined in Eq. (6) as G(t) = (ra(T)a^(O)), is plotted in Fig. 10 for two set of values of 
the parameters {J/U^ f^/U), one in the MI phase, the other in the SF, both very close to the tip of the first insulating 
lobe. 

Let us first remark that the decay of the Green function is independent on temperature in this temperature range, 
the effect of temperature being only to cut the decay around |r] ~ (3/2 due to periodicity. Therefore we can safely 
assume that the Green functions reported in Fig. 10 are good estimates of their zero-temperature limit. 
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FIG. 9: Interaction dependence of the transition temperature at unit filling. The result obtained on the c = 6 Bethe lattice 
(black lines and filled circles) is compared to Monte Carlo simulation (blue lines and open diamonds) on the cubic lattice (data 
from [18]), and to the prediction of the mean-field approximation (red dashed curves). 

For the interpretation of these curves it is worth recalling the spectral representation of the Green function [35] at 
zero temperature: for r > 0, G{t) = \ (a) [^ + /q pp{E) exp{—ET), where Pp{E) is the (on-site) spectral density for 
particle excitations. A similar expression holds for r < in terms of hole excitations. The long time limit of G{t) is 
thus strictly positive in the SF phase, and we checked that is in agreement with the value [ (a) [^ previously computed 
(see the dotted line in the right panel of Fig. 10). 

Let us now turn to the discussion of the Green function in the MI. In the mean-field description of this phase the 
hopping is completely suppressed, hence particle (resp. hole) excitations corresponds to the addition (resp. removal) 
of one particle from the background of the commensurate filling of immobile particles. To be concrete let us consider 
the first lobe with (rt) — 1; the energy of these particle and hole excitations are easily found to he Ep = U — p and 
Ef, = ^, hence follows the expression of the Green function [45]: 

G„,f (r) = 0(T)2e-(^-^)" + e{-T)e-^^ . (35) 

In other words the spectral density of particle excitations is made of a single delta function (and similarly for the hole 
excitations). The mean- field Green function is reported in Fig. 10 as a dashed black line; the comparison with the 
results of the Bethe lattice computations shows that the latter is more complex. The decay of the latter is not a simple 
exponential, thus reflecting a non-trivial spectrum of excitations. Moreover the energy of particle and hole excitations 
is lower on the Bethe lattice, as can be seen from their slower decay. This happens because on the Bethe lattice, due 
to the local finite connectivity, particles can still hop around even in the MI phase, resulting in a gain of kinetic energy 
and a lowering of the energy cost of the excitations. Note that different improvements of the mean-field treatment, 
for instance Random Phase Approximations (RPA) [13, 14], provide a much better description of the spectral density 
in the Mott phase, in particular a good approximation to the wave-vector dependence of the excitation energies for 
finite-dimensional lattices. A quantitative comparison with the predictions of the RPA is beyond the scope of this 
paper. We focused instead on the comparison with the fully-connected result and B-DMFT, because they are related 
to the Bethe lattice via a well-controlled limit, namely the large connectivity one. 

We did not attempt to perform the inverse Laplace transform to determine the spectral function, yet we performed 
the following analysis to estimate the scale Ep of the slowest decay in the r > part of the Green function. A plot 
of log[r"G'(r)]/r is seen to approach a constant at large r for a '--^ 1, suggesting that G(t) decays as exp(— i?pr)/r, 
hence that the spectral function Pp{E) is finite a.t E = Ep and vanishes for E < Ep. This is consistent with the 
results of [13, 14]. For simplicity, we have chosen to fit the Green function to G(t) = 2e^^^'^{l — e~^p '^)/(Ap r), 
corresponding to a flat density of states in [Ep, Ep + Ap]. We stress that this is an arbitrary choice that we made only 
to fit the long time behavior of G(t) and determine Ep. Note that at very small r <C 1/J, hopping is irrelevant and 
G(r) ^ Gnif (t) at first order in t. This imposes a relation between Ap and Ep, hence there is a single free parameter 
in our fitting procedure. 

In Fig. 11 we show the behavior of the Green function on approaching the phase transition from the MI phase, 
close to the tip of the first lobe by varying J/U at constant p,/U = 0.39. We focus on positive r and fix /?J = 6, 
which according to Fig. 10 is a low enough temperature, such that the Green function is identical to its zero- 
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FIG. 10: Green function G(r) as a function of the imaginary time at c = 4. (Left panel) Results in the MI phase, at 
J/U = 0.0523, fi/U = 0.39, A'^cut = 4 and l3J = 2 (black dashed-dotted curve), 4 (red dashed curve), and 6 (blue continuous 
curve). Also reported is the mean-field result for the same parameters (black dotted curve). (Right panel) Results in the SF 
phase, at J/U — 0.0555, n/U = 0.39, A'cut = 6 and /3J — 2 (black dashed curve) and 4 (blue curve). Also reported is the long 
time limit of G(r), which corresponds to the value of (a)^ for these values of the parameters. 



n 





1 


2 


3 


4 


5 


P„ (SF) 


0.034 


0.93 


0.034 


0.00027 


5.7e-07 


2.2e-ll 


Pn (MI) 


0.023 


0.95 


0.023 


0.00013 







TABLE I: Values of P„ in the MI and in the SF, for the same values of the parameters as in Fig. 10, namely J/U — 0.0523, 
fi/U = 0.39, iVcut = 4 and /3J = 4 (in the MI phase) and J/U = 0.0555, fi/U = 0.39, iVcut = 6 and /3J = 4 (in the SF phase). 



temperature limit for r < /3/2. We use the rescaled time variable tU, in such a way that the mean field result 
G'nif(''") = 2exp[— rJ7(l — fi/U)] yields the same decay at all values of J/U. On the contrary, the Bethe lattice Green 
function depends quite strongly on J/U, and its decay slows down on approaching the transition. The result of the 
fitting procedure explained above is shown in the inset of Fig. 11 and shows that the parameter Ep is smaller by a 
factor of 2 on the Bethe lattice compared to the mean field results; it decreases on approaching the transition, but 
remains finite at the transition. Indeed we have computed here the on-site Green-function, which reflects only the 
localized single-particle excitations, whereas the transition towards the superfluid phase is towards a delocalized state 
(small momentum in the finite dimensional case [18]). To detect the signature of the transition in such a way it would 
be necessary to compute the spatial and temporal dependence of the Green function, which is in principle possible 
using the cavity method (see next subsection). 

We expect (but did not check in detail) that also in the SF phase the on-site Green function will decay exponentially 
up to the transition point. This shows also that the BEC state on the Bethe lattice has a peculiar nature as compared to 
the finite dimensional case: condensation happens in the uniform state, which is separated from the other eigenvalues of 
the connectivity matrix by a gap [46, 47]. In particular, one can check explicitly by a Bogoliubov-like computation [48] 
that the single-particle excitation spectrum is gapped at small U/J. 

The difference between the mean field approach and the Bethe lattice result is also unveiled by the analysis of the 
probability distribution P„, already studied recently in [19] using Monte Carlo simulations, defined as the probability 
to detect n bosons on a given lattice site: 

Pno = (^«,no> = |Tr [5„,„„e-''^] . (36) 

This is another quantity sensitive to the occupation number fluctuations. In table I we report the values of P„ in the 
MI and in the SF, for the same values of the parameters as for Fig. 10. These data clearly indicate that even in the 
MI phase the occupation number of a given lattice site is not strictly fixed to an integer value, as in the mean field 
approach: there is still a finite probability of finding a number of bosons different from one within the lobe at (n) = 1. 
Interestingly enough, these quantities can be probed in experiments [49]. 
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FIG. 11: (Main Panel) Green function G{t) as a function of the rescaled imaginary time tU for c = 4 in the MI phase, close 
to the tip of the first lobe at n/U — 0.39, A^'cut ~ 4 and PJ — 6, at different values of J/U (approaching the transition to the 
superfluid, from bottom to top). Also reported is the mean-field result for the same parameters (black dashed curve). (Inset) 
Particle excitation gap Ep/U is plotted as a function of J/U (see text for details); the dashed black line is the mean field result 
Ep/U = 1 — IJ./U, the vertical line marks the transition to the superfluid. 




4- One particle density matrix 

We present now our results for the one-particle density matrix pij = ^ajflj^ Because of the invariance structure 

of the Bethe lattice this will depend (in the thermodynamic limit) only on the distance between the sites i and j, 
that we denoted as d{i,j) in the discussion of Sec. Ill A. In Fig. 12 we report at different points in the {J/U, n/U) 
plane, crossing the insulator-superfluid transition around the tip of the first lobe. 

In the MI phase pd decays to zero, while in the SF phase we checked that it correctly decays to | (a) p. The decay is 
always quite fast, and in the SF phase we do not have enough points to perform a reliable fit to extract a correlation 
length, altought it is clear that the decay becomes slower and slower on approaching the transition to the MI. On the 
contrary, in the MI phase it is very easy to extract the correlation length by measuring the slope of log pd versus d 
for d > 3. The result is plotted in the inset of Fig. 12 and shows that the correlation length increases also on the MI 
side. 
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It might be surprising at first sight that the correlation length stays finite at the transition; however, as we discussed 
at the end of section III A, this is a general property of Bethe lattices that is related to the fact that the number of 
sites Md at distance d from a given one scales as (c — l)'' for large d. Despite this difference at criticality, we believe 
that the possibility of investigating spatial correlations on the Bethe lattice is very interesting away from critical 
points, and is an important improvement of the cavity method with respect to the mean field approximation. 



V. DETAILED DERIVATION 



A. Lattice boson models in the Suzuki-Trotter formalism 



This section is devoted to a complete derivation of the equations used in Sec. IV. We shall consider a slightly more 
general Hamiltonian than (1), 



N 



H = 



(37) 



i=l 



where the hopping strength J{ij) is allowed to vary from edge to edge and we consider arbitrary local potential 

energies Vi{n). This will be convenient both for notational reasons and because part of the following discussion will 
also apply to the disordered Bose-Hubbard model. The original model (1) is recovered by taking J{ij) = J on all 
links, and Vi{n) = ^n(n — 1) — nn. The Hilbert space of the model is spanned by the occupation number vectors 
. . . , njv), where is a positive or null integer counting the number of particles on site i. We shall use the more 
compact notation \n) = \ni, . . . ,nAr) to denote one of these basis vectors. For completeness we recall the action of 
the annihilation a, and creation at operators in this basis. 



(38) 



ai\n} = ^/r^\nl, . . . ,ni-i,ni - l,ni+i, . . . ,nN) , 
al\n) = ^/r^ + l\nl,...,ni--l,ni + l,ni+l,...,nN) ■ 



The number operators aja^ arc diagonal in this basis, with ajaijn) = ni\n). The partition function of the Hamiltonian 
(37) will be computed using the Suzuki- Trotter formula (18), cutting the imaginary time axis of length /3 in a number 
Ng (ultimately sent to infinity) of slices. At each of these slices a representation of the identity in the occupation 
number basis is inserted. This yields 



Z = Ti[e 



-f^"] = lim 



exp 



N 



a=l i=l 



a=l 



(39) 



The index a is the discrete coordinate between 1 and A^s along the imaginary time axis, and we use periodic boundary 
conditions n^^"*"^ = n}. The quantum problem has thus been transformed in a classical one, in terms of the imaginary- 
time dependent classical variables n". The expression above of the hopping interaction is however unpractical for our 
future needs, we shall thus transform it by introducing a set of auxiliary variables = {yfi^j} which can take values 
or 1, depend on the discrete time a, and are defined on the directed edges of the graph (an edge thus bears, 
for each time a, two variables yf_^j and Uj^i)- It is simple to show that 



■ <i.J) 



En 

y" {hi) 



(40) 



where here and in the following 1(A) = 1 if the condition A is fulfilled, otherwise. We take by convention a;^^° = 1 
for any value of x (including zero). The justification of (40) can be done by inspecting the behaviour of its left and 

right-hand-side order by order in l/Ag. The leading term corresponds to = n""*"""^, and indeed all y's must vanish 
at this order. At order l/Ag a single y, say yiUj, is equal to 1, moaning that a boson has hopped from site i to 

site j between the discrete times a and a + 1. The term \ n^^^nf follows from the action of creation/annihilation 
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operators recalled in (38). The kinetic energy part of the partition function (39) can be transformed by using the 
representation (40) for each time a; reordering the various terms leads to a compact expression, 

N 

Z = ^lini^ n ^(^J> (y(*,i>) n {y{i,j)}jedi) ■ (41) 

Following the convention we already used in Sec. Ill bold symbols stand for time dependent quantities, with a discrete 
time coordinate a e [l,iVs] as long as A^s is finite, or a continuous time r £ [0,/3] in the limit Ng oo. The 
correspondence between the two notations is r = a[3/Ns. We also use the shorthand y". --^ = {yf_^j,yj'^i) for the pair 
of hopping variables on both directions of an edge The weight of a configuration (n, y) factorizes into a product 

over all edges of 

-(M)(y(M))=n(^) ' (42) 

a=l ^ ^ ^ 

the hopping strength J{ij) is thus conjugated to the number of hopping events between sites i and j. The second 
part of the weight in Eq. (41) is a product over the sites of 



Wiirii, {y{ij)}jedi) = exp 




The first line arises directly from the local potential energy term of the Hamiltonian. The second line enforces the 
consistency of the occupation number trajectory nf with the hopping events that occur on the adjacent edges, and 
incorporates the ^/n terms arising from the action of the creation/annihilation operators. 

The expression of the partition function in terms of a classical model with imaginary time-dependent variables we 
just obtained is valid for any underlying graph. We want to emphasize that the spatial structure of the classical model 
is the same as the original one: the variables rii are located on the vertices of the graph, the y/^^ on the edges, and 
in the weights Wi the interactions only occur between edges and site variables which were adjacent in the original 
model. In particular, if the latter was locally tree-like, this remains true for the representation (41) on which the 
Bethe approximation can be performed, as we shall do in the next subsection. 



B. The Bethe free energy 

One can follow two roads to obtain the self-consistent equation (27). The first one corresponds to the reasoning 
presented in the simpler case of the Ising model in Sec. Ill A. Assuming the graph to be a tree, one can easily compute 
recursively Zi^j{y(^^j^), the partial partition function for the subtree rooted at i excluding j, for a given value of the 
hopping trajectory y(^ijy Introducing its normalized counterpart r}i^j{y ^^ ^s^) = Zi^j{y ^^ ■'^)/ Zi^j, one obtains 

^i^i{V{i,3)) = ^.^{i,3){y{i,3)) H w^i("i>y(i,j>'{y(i,fc>}feeaAj) n ^kMy{i,k)) ■ (44) 
The probability distribution of the occupation number trajectory on site i is then expressed as 

Vi{ni) = ^ Wi{ni,{y{i,j)}jedi)]Jvj^i{y{i,j)) ■ (45) 

* {l/<i,j>}je8i oedi 

These last two equations are the analogs of Eq. (13) of the Ising model. If one considers now the homogeneous case 
of a Bethe lattice of connectivity c, with the same hopping strength J{ij) = J on all edges and Vi{n) — V{n) for 
all sites, the probability laws r]i^j take a common value Tjcav on all directed edges, which is the solution of Eq. (27). 
In that equation y is an hopping trajectory on an edge, characterized more explicitly as y = (j/^,yl, . . . jU^" ,y_''), 
where = 1 (resp. = 1) means a particle has arrived (resp. leaved) the considered site at discrete time a. 
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The expression of the weights wunk and wnei of Eq. (27) is found by specializing (42) and (43) to the homogeneous 
situation, 



(3J 



e a=i 



wun^iy) = ^^J-' , (46) 

and 

waiter (y,yi---,yc-i) = X!' 

n 

n |l (n"+^=n"+y'^-yl + J2[yh-y-,i]^ (v^)''"^^''^' | . (47) 

Note that the direction of y^, for i € [1, c — 1], is reversed with respect to y, hence the inversion of the signs on the 
variation of n" from one discrete time to the next. Finally the weight Wgite used in Eq. (32) to obtain the probability 
distribution of an occupation number trajectory reads 

ti^site(n,yi...,yj = e " n |l U"+i = n« + ^^[y?., - j (v^)--'^'" (v^^ (48) 

The second strategy to reach the equation (27) on r/cav consists in first writing the Bethe approximation [39] for 
the free-energy of the model defined in (41) for an arbitrary graph. This yields 

11^/ 
F = --lnZ= - -^In ^iini,{y{i,j)}jedi)'[[^J^iiy{i,j)) 

i=i \n,,{y(ij}}j€di 

One can check that this relation is exact whenever the graph is a tree, otherwise it corresponds to the Bethe ap- 
proximation. An important property of this expression of the free-energy is its variational character: the stationarity 
conditions with respect to the parameters rji^j are nothing but the local recurrence equations (44). This means that 
the computation of the derivatives of the free-energy with respect to its parameters /x, (3, J{ij) , ■ • ■ can be performed 
by considering only the explicit dependence on these parameters, discarding the implicit dependence through the 
rii^j. We shall also show in Appendix A how to devise a further (static) approximation based on this variational 
expression of the free-energy. Considering the particular case of an homogeneous Bethe lattice leads to the following 
free-energy per site: 



/= lim -J-lnZ-- 

■' N^oo /3N 13 



^ In ^ E ^site(n, yi, . . . , Vc)'ncAVl) ■ ■ ■ »?cav(yc) j + ^ ^^^^^(^) ^cav(y) 



(50) 

where we have used the fact that a regular graph of N vertices and connectivity c has ciV/2 edges. One can see that 
(27) is the stationarity condition of (50) with respect to r^cav 

By using the recurrence equations (44) it is possible to rewrite the Bethe approximation of the free-energy, for an 
arbitrary graph, as 

N 

F = -1 InZ = -1 ^ In i^Zi^^z^^i) + i ^ ^ In^, , (51) 

where Cj is the degree of vertex z, and the various z^s correspond to the normalizations in (44), (45). This form does not 
have the variational property explained above; it has however the advantage of being easier to compute numerically. 
In particular for the homogeneous case it yields 

where z^av and Zgite are the normalizations in, respectively, Eqs. (27) and (32). 



21 



C. The resolution of the cavity equation 

As we have explained in Sec. IV A we use a weighted sample representation of r/cav to solve (27). The point we 
want to specify more precisely here is the method used to sample hopping trajectories from the conditional law 
-P(y|yi) • • • ) Vc-i) ^^'^ to compute the associated normalization -Z(yi, . . . , Vc-i) (both defined in Eq. (31)). A similar 
construction can be found for the slightly simpler case of quantum spin 1/2 models in [34]. 

Let us begin with the computation of Z. A first important remark, already mentioned above, is that in the 
continuous time limit {Ng oo) the hopping trajectories y typically contain only a finite number (with respect to 
Ns) of hopping events, i.e. of discrete times a where ^ 0. We can thus assume without loss of generality that 
the hopping events for the different trajectories occur at different values of the discrete time. Let us call p the total 
number of hopping events occuring in (y^, . . . , j/c-i): denote ai < ■ ■ ■ < ap their discrete time of occurence. We 
consider the Hilbert space of a single site, with a and the annihilation/creation operators, and call bj = a (resp. 
bj = a^) if the hopping of time aj is towards (resp. outside) the vertex under consideration. Finally we define = bj 
when a = aj, Ca = 1 (the identity operator) otherwise. A moment of thought reveals that 

2{yi, ■ ■ • ,yc-i) = E^'-k(y)«;iter(2/,yi, • • • ,ye-i) = ^ n<""l^^'"'''"*"^^"'^"^"*"^«l""^') > (^3) 

y n a=l 

up to corrections of order Ng^'^ . The sum over n can thus be written as a trace of this product of operators. To perform 
the continuous time limit it is convenient to define Tj = -^oij, which are the continuous times of the hopping events 

in (y^i • • • ,yc-i)- We also introduce W(X) = e-^(-^("^'')+-^(<»+''^)), the propagator of an imaginary time evolution on 
an interval of length A for a single site Hamiltonian V{a)a) — J{a + a^). We thus obtain finally 

Z{y^,. . . ,y,_i) = Tr (w{n)biW{T2 - n)b2 . ..W{Tp - Tp.i)bpW{P - r^)) . (54) 

The single site Hilbert space is a priori infinite dimensional because the number of particles is not bomided. Very 
high occupation numbers are however unlikely because of their prohibitive potential energy. We can thus safely put 
a cutoff A^cut on the possible values of n, which will not change the properties of the model if A^cut is sufficiently large 
with respect to the average density of particles. This would amount formally to take V{n) = +oo for n > A'cut- In 
this case the dimension of the Hilbert space is reduced to A^cutj and (54) is nothing but the trace of the product of a 
finite number of finite size matrices, whose numerical evaluation does not present any particular difficulty. 

We turn now to the problem of the generation of an hopping trajectory y given the ones of the c — 1 other 
neighbors, according to the law P{y\yi, ■ ■ ■ , yc-i) stated in (31). As explained in Sec. IV A (see in particular Fig. 5), 
we can determine y by first drawing the occupation number trajectory n and then deduce y from its discontinuities 
not associated to hopping events in {y^, . . . , i/c-i)- We stick to the notation n < • • • < Tp and bi,. . . ,bp for the 
parametrization of the hopping events in (j/^, . . . , y^-i). Let us call no = n(r = 0), (resp. n-) the value of n(T) at a 
time just after (resp. just before Tj+i), with the conventions n'p = hq. The joint probability law of these occupation 
numbers which arise from the expressions of wunk and Witei given in Eqs. (46), (47) reads in the continuous time limit 

P{no,n'o,ni,n[, . . . ,np\yi, . . . ,y^_i) = — -(no|W^(Ti)|no) TT |(n-_i|6i|ni)(ni|W^(Ti+i -Ti)|n-)| , (55) 

^lyi) • • • )yc-i; 

with Tp-i-i = (3. This probability law is well normalized according to the above expression of Z{yi, . . . , y^._i). Moreover 
because of the "imidimensional" structure of the imaginary time axis it is rather simple to generate a set of integer 
numbers (no, nf,, ni, n'^, . . . , Up) according to this law. Once these intermediate occupation numbers are known, the 
generation of the occupation trajectory n(r) can be done independently on each of the p + 1 time intervals between the 
imposed hopping events. On the z+ I'th interval between and Ti+i one has to generate the trajectory corresponding 
to a path integral representation of the effective Hamiltonian V{a^a) — J{a + a^), conditioned to begin in rn at time 
Ti and to end in n[ at Ti^i. 

To simplify the notation let us consider this problem for an interval of imaginary time [0, A], with boundary 
conditions n(0) = n, n(A) = n'. One way to justify the sampling procedure is to start with the identity 



JO 



(56) 



22 



valid for any non-commutative operators X and Y. We apply it with X = —V{a^a) and Y = J{a + a^), which yields 



A 



{n\W{X)\n') = (n|e-^^("^") + / dr e'^^^'^^'^) J(a + a^)t?(A - T)|n') (57) 

/■A _ 

= 5n,n'e-^^^"^ + J / dr e-^^(")(n|(a + a^)W{\ - r)|n') . (58) 
Jo 

The interpretation of this formula is as follows: a path representative of the matrix element (n|M^(A)|n') is either a 
constant trajectory (possible only if n = n'), or it is constant up to a time r, where it jumps to n ± 1, and is followed 
by a path of length A — r with boundary conditions n± 1 and n'. We can thus sample the path n(r) according to the 
following recursive procedure: 

• if n = n' and with probability e~^^^"'^ /{n\W{X)\n), set n(r) = n for r e [0, /3] and exit the procedure 

• otherwise 

— draw a time r e [0, A] with the probability law written below in Eq. (59) 

— draw a = 1 with probability ^ , a = -1 otherwise 

— set n(r') = n for r' e [0,t] 

— call the same procedure with initial condition n + a, final condition n', for the imaginary time interval [r, A] 

The probability law for the time t of the first discontimiity in the number occupation trajectory follows from (58); it 
is more convenient to express it as the probability G{u) that t < u, 

/;dT e-^(")(n|(a + at)T?(A-T)K) 
G(u) = -^3; ^i^^ . (59) 

J^dT e-^nn){n\{a + a^)W{X-T)\n') 

Because of the cutoff on the possible number occupation one can numerically compute this function by diagonalizing 
the effective Hamiltonian (only once for each choice of the parameters U,fi,J). Drawing r then amounts to use a 
random number G uniformly in [0, 1] and set r = u~^{G) (the functional inverse of G{u)). 



D. The computation of the observables 

We will show now that, as announced in Sec. IV A, all the physical observables can be obtained as averages over 
r/cav of some functions of the hopping trajectories. The computation of such an average can then be done as a simple 
sampling of the population representing ?7cavj as explained in Eq. (34). 

Let us begin with the free-energy. Using the formula (52) its computation amounts to the determination of the 
normalization constants Zcsev and ^site- The former reads 

^cav= XI ^cav(yi)---?7cav(yc-l)2^(yi:---.yc-l) - (60) 

l/l.->l/c-l 

which is readily evaluated since we obtained an explicit form of 2{yi, . . . ,yc-i) ^Q- (54)- The latter is seen from 
(32) to be 

Zsite= X] '?cav(2/l)---%av(2/c)X]^^'*«*("'yi'---'^c) • (61) 

Using again the notation ti < • • • < Tp and bi, . . . ,bp for the paramctrization of the hopping events in {yi, . . . , y^) 
and exploiting the expression (48) for Wsite, we obtain in the continuous-time limit 

^site= '?cav(yi)... r?cav(yc)TV (H^o(ri)6i W^o(r2 - ri)62... 6pWb(/3 - Tp)) , (62) 

where we defined Wo(A) = g-^^f" '^), the single-site propagator without its hopping term. Indeed all the possible 
hopping events to and from the considered site are fixed by the c trajectories y^, . . . , y^. This trace is over a product 
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of finite size matrices (thanks to the bound A'cut on the occupation numbers) and hence computationaUy harmless. 
Note that only configurations of y^^ . . . which encode the same number of jumps towards/outside the central site 
do contribute, because Wq is diagonal in the number basis. 

We consider now the thermal average ^g'(ajai)^ of an arbitrary function q{n) of the number operator on one site. 
Prom the expression (32) of the probability distribution r]{n) we obtain 

q(alai)) = — V rfcav(yi) ■ • ■ r?cav(yc)Tr f g(aMW^o(Ti)6iVKo(r2 - ti)62 ■ • ■ W(/3 - Tp)) . (63) 

This formula allows to compute the density of particles, the average local energy ei, and the occupation probability 
Pno by taking q(ri) = n, q{n) = V{n) and q{n) = Sn,noj respectively. 

The order parameter (aj) is obtained by the insertion of an annihilation operator in the effective single-site problem, 

{ai) = — V 7?cav(yi) . • . %av(2/c)T> (aWb(Tl)6l Wo(t2 - Ti)62 • • • bpWo{P - Tp)) . (64) 

One has indeed by definition {ai) = Ti:[aie~^^]/Z , where the trace is here over the Hilbert space of the A^-sites 
Hamiltonian. One can reproduce all the steps leading to the representation (41) of the partition function, with the 
additional operator Oj modifying the expression of the weight Wi on the considered site. The other sites are however 
left unmodified, and hence integrating over them leads to the same equation on Tjcav The presence of aj only shows 
up when the sum over the degrees of freedom of site i is performed, and leads to this insertion of the annihilation 
operator in the single-site computation (64). We now come back on the problem of the initial condition in the 
population dynamics mentioned in Sec. IV A. One can indeed see that (64) strictly vanishes whenever all hopping 
trajectories in the support of Tjcav have the same number of jumps in the two directions of their edge. Moreover this 
symmetry is conserved by the iterative equation (27); we thus had to initialize the population dynamics algorithm 
including asymmetric hopping trajectories. In the "insulating" phase these trajectories disappear during the iterations, 
while in the "superfluid" a finite fraction of them keeps the asymmetry, thus allowing for a non zero value of (ai). 

The computation of the Green function G'>(r) = Tr[e~(^~'^)^aie^'^^aJ]/Z can be done similarly, with now the 
insertion of two creation/annihilation operators in the single-site problem. Similarly to Eq. (63), and using the same 
notations, one obtains 



. {r) = — ^ ??cav (y 1 ) • • ■ ?7cav (Vc) > 

X Tr (aWo(Ti)6i Wb(r2 - n)b2 • • . hWoir - n)a^Wo{n+i - T)bi+i ■ ■ ■ bpWo{l3 - Tp)) , 



-^i.-.i/c (65) 



where the index i is such that Tj < r < Tj+i. 

We turn now to the computation of the average kinetic energy ex = {Hk) /N. By definition it is equal to the 
derivative of the free-energy per site with respect to J. Using the variational property of the expression (50), we can 
thus write 

E^;;ii(io'7cav(y)|y| 

V 

where we have defined \y\ as the total number of hopping events in y. Indeed the only explicit dependence on J in 
(50) is in the weight wunk given in Eq. (46). Replacing one of the Tjcav by its expression (27) both in the numerator 
and the denominator, and reordering terms leads to 

eK = -^^ ^cav(yi)...??cav(yc)|yi|Tr(^?^o(Ti)6iWo(r2-Ti)62...6pW^o(/3-Tp)^ . (67) 

Finally we consider the computation of the one-particle density matrix pd for two sites at distance d. For notational 

simplicity we shall call them and d, pd = (^aljCid^- The unique path linking these two sites is schematized in Fig. 13; 

we have distinguished the hopping trajectories y^ along the chain for i e [0,d — 1], from the contribution of the 
edges outside the chain, y- The index j runs from 1 to Cj, with cq = ca = c — 1 and ci = • • • = Cd-i = c — 2, 
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FIG. 13: Computation of the one-particle density matrix at distance d, pd = \ ajad 



as the extremities and d have one less neighbor in the chain. The computation of ^aja^y is a generahzation of 

the determination of (oj) explained above: one writes ^aja^^ = TT:[a\ade~^^\/ Z and expresses the numerator with 

the Suzuki- Trotter formula where the local weights Wq (resp. Wd are modified with the insertion of a creation (resp. 
annihilation) operator. Summing over all degrees of freedom outside the chain of sites between and d leads to 

%ad^ E n»?cav(y-j) E n^i'iink(yi) E Wsite{no,yody'oj})Uwsite{ni,yi_^,yi,{y'^j})wsite{nd,yd-^^ ' 

(68) 

where in the numerator are modifications of the vertex weight (48) due to the insertion of the creation/annihilation 

operators on sites and d. Integrating progressively the degrees of freedom along the chain one can easily shows 
that the denominator of (68) equals z^^^Zsitc- To compute the numerator we shall define a sequence of probability 
distributions ?7cav, as the law of the hopping trajectories along a chain where the part between i + 1 and d is 
removed. This can be computed by recurrence on i, with 



^cal(y) = 7^^^1ink(y) J2 ^cav(yi)---»?cav(yc-l)X!^ite(".y'yi.---.yc-l) > (69) 

and 



(0) 

l/l>->l/c-l 



-^cav 



^Ll^^y) = ;^(iq:Ty^link(y) J2 ^clUyi)^cav(y2) ■ • •»?cav(yc-l) X^^^'t^^^'^'yi' • • • '^c-l) • (70) 

i/i>---,yc-i 



-^cav 



These probability distributions can be encoded numerically as weighted samples, in the same way as we did for jjcav 
Once they have been determined up to distance d — 1 we can compute 

4fte = Yl '^cf^^^ (yi)»?cav(y2) • • • ??cav(yc) X] ^s7te(^. yi> • • • > yc) • (71) 



One then finds 



Jd) d-1 U) 
'^c-wr^ T r 



VI. CONCLUSIONS 



In this paper we have extended the analytical tools of the quantum cavity method to Bethe lattice bosonic models. 
The cavity approach yields a self-consistent equation for the probability of hopping trajectories in continuous imaginary 
time, which can be efficiently solved numerically using the population dynamics algorithm. Though this numerical step 
is unavoidable to obtain quantitative predictions it is distinct from a Quantum Monte Carlo algorithm, in particular 
the thermodynamic limit is taken analytically. 

We have shown how this method allows to compute easily several physical observables such as thermodynamic 
quantities (average occupation number, energy, free energy, condensate fraction, . . . ) as well as imaginary time Green 
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functions. Let us emphasize in particular that computing the free-energy (which is often difficult with other methods, 
requiring for instance thermodynamic integration in Monte Carlo simulations) is crucial to determine the location of 
first-order phase transitions which can occur in generalized versions of the model. 

The finite connectivity of the Bethe lattice leads to two kind of improvements with respect to the mean-field 
description: a distance between two sites can be defined, which allows the computation of spatial correlation functions. 
Moreover the description of the Mott Insulator phase is richer, with a non-trivial effect of the hopping of the particles. 
The prediction for the zero-temperature phase diagram of Bethe lattices with connectivity 4 and 6 is in reasonable 
agreement with the Quantum Monte Carlo simulations in 2 and 3 dimensions [17, 18]. 

Let us sketch some possible directions for future work. A first possibility would be to explore the next levels 
of the Cluster Variation Method [26] hierarchy of approximations, which contains the mean-field and the Bethe 
approximations as its first two steps. This approach, similar in spirit to the Cluster DMFT [50], should improve the 
accuracy of the predictions for small dimensions. It should also be interesting to investigate the effect of disorder in 
the Bethe lattice model; if the Bose Glass phase [1, 27] is destroyed at the mean-field level, it could be possible to 
describe it on the Bethe lattice, which can indeed display localization properties [28]. Another open direction is the 
study of several generalizations of the Bose-Hubbard model, including bosons with spin, multiple interacting species 
of particles [23, 24, 51, 52], and interactions between nearest- neighbor sites [53]. In particular in the latter case, 
we expect that the presence of interactions inducing geometrical frustration will give rise to glass phases, since the 
same happens in the classical limit of zero hopping [54] . It would be very interesting to check whether these glass 
phases could exhibit Bose-Einstein Condensation as it happens for their crystalline counterparts [55] . This would add 
some new element to the ongoing discussion on the relevance of disorder to induce superfluidity [56], and might be 
important to interpret and understand recent experiments on supersolidity [57]. 

We finally point out that, although in this paper we have only focused on the very low temperature regime (except 
for the transition temperature at unit filling displayed in Fig. 9), the method presented here can be applied at higher 
temperature (and is even numerically easier in that case), which can be useful for the interpretation of cold atoms 
experiments [58]. From this perspective the case of unidimensional systems should deserve a special attention, as 
they fall in the category of tree structure and could thus be solved with the cavity method, even in the presence of a 
site-dependent chemical potential modelizing the trapping field. 
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APPENDIX A: STATIC APPROXIMATION 



We have obtained in (49) a variational expression of the Bethe free-energy, in terms of the parameters rii^j{y j^). 
A further approximation consists in restricting the rji^j to a simply parametrized subspace, and finding the extremum 
of (49) within this simpler ansatz. A possible choice for this parametrization, inspired by the so-called static approxi- 
mation developed for quantum spin models [34, 59], consists in retaining in Vu^j) only the information on the number 
and direction of the jumps, but assuming their times of occurence to be uniformly distributed. This corresponds in 
formula to the ansatz 



(Al) 



where pi^j is the distribution probability for the number of jumps from % to j and viceversa. The two terms of the 
variational free-energy (49) can be computed within this ansatz. One obtains indeed in the limit — > oo: 

1 °° ! ! 
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and 

(A3) 

where 



Ai{N+,N_) = 



-Tr 



-/3(Vi(a^a) — z+a^ — 2_a) 



(A4) 



z+=z—=0 



As a simple illustration of this static approximation one can recover the mean-field free-energy (9) for the case of a 
Bethe lattice of connectivity c with uniform hopping strength J and the same potential energy V{n) on all sites. This 
is indeed obtained after a short computation by taking a Poisson distribution of average /3JA for the number of jumps 
in each direction, 



p{n+,n-) = e 



-2g7A (/?JA)"++"- 



n+!n_! 



(A5) 



where A is proportional to the parameter tjj in (9). 
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The initial condition should however allow for the possible spontaneous symmetry breaking (a) 7^ 0; we shall discuss this 
point further in Sec. V. 

Note that we compare our results for a Bethe lattice with finite connectivity with those obtained in [24] using the semi- 
circular density of states. The latter is strictly valid in the limit of large connectivity, but this effect should be of subleading 
order in the large connectivity expansion. 

In particular the cusp in the kinetic energy at transition from the SF to the MI becomes sharper and ex becomes slightly 
smaller as we go from 0J = 1 to /3J = 2 in the first superfluid region. No further modifications are observed within our 
numerical errors for f3J > 2. 



